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ABSTRACT 

This paper describes SUNRISE, a parallel, free Monte-Carlo code for the calcula- 
tion of radiation transfer through astronomical dust. SUNRISE uses an adaptive- 
mesh refinement grid to describe arbitrary geometries of emitting and absorb- 
ing/scattering media, with spatial dynamical range exceeding 10 4 , and it can effi- 
ciently generate images of the emerging radiation at arbitrary points in space. In 
addition to the monochromatic radiative transfer typically used by Monte-Carlo 
codes, sunrise is capable of propagating a range of wavelengths simultaneously. 
This "polychromatic" algorithm gives significant improvements in efficiency and 
accuracy when spectral features are calculated, sunrise is used to study the 
effects of dust in hydrodynamic simulations of interacting galaxies, and the pro- 
cedure for this is described. The code is tested against previously published 
results. 

Key words: dust - radiative transfer - methods: numerical. 



1 INTRODUCTION 

Interstellar dust profoundly affects our view of the uni- 
verse, from obscuring the stars forming in giant molecular 
clouds in our Galaxy, to camouflaging extreme starbursts 
as relatively unremarkable galaxies in Ultraluminous In- 
frared Galaxies, to completely hiding high-redshift galax- 
ies from our view except in the infrared, as in the sources 
detected with SCUBA. Because of this, modeling the ef- 
fects of dust has been the subject of an ever-increasing 
number of papers. Initial models used very simple as- 
sumptions, such as the dust being distributed as a fore- 
ground screen. While appropriate for observations of sin- 
gle stars, this assumption fails miserably in the case of 
galaxies, where the stars and dust are intermixed. In this 
situation, the appearance of the system will depend not 
only on the characteristics of the dust itself but also on 
the relative distributions of stars and dust. In this sce- 
nario analytic solutions in general do not exist, and nu- 
merical solutions become necessary. 

Numerical approaches to solving the radiative- 
transfer problem for dust can gener ally be classified 
as either finite-differen ce methods (e.g. ISteinacker et alJ 
l2003l : lFolini et al.ll2 003'l . or Monte-Carlo methods, where 
the problem is solved in a stochastic sense. The advan- 
tages of Monte-Carlo methods are that they can easily 
treat complications such as arbitrary geometries and non- 
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isotropic scattering, the main drawback being that gen- 
erating a sufficiently accurate solution can be computa- 
tionally quite expensive. Traditionally, the computational 
expense has been reduced by using simple geometries and 
exploiting symmetries in the problem, such as assum- 
ing spherical or azimuthal symmetry. It is only more re- 
cently that fully th ree-dimensional, arbitrary geometries 
have been tackled fWolf et alJll99St iBianchi et al.ll200Cl; 



Gordon et alJl2001l:lKurosawa fc HillierlbOOltlBaes et alJ 
20031 : lHarries et al.l [2004) . Codes combining dust and 

ffiotoionization radi ative transfer are also appearing 

Ercolano et al.ll2005l) . 



This work describes sunrise, a new code for Monte- 
Carlo radiative transfer calculations. The development 
of this code was motivated by a desire to calculate 
the effects of dust in a large program of hydrody namic 
N-body simula tions of interacting galaxies (ICox 2004; 
ICox et alJl2005l) . The main requirements for this are that 
the code be able to handle arbitrary geometries and 
resolve the large spatial dynamic range in the simula- 
tions, and that the large number of calculations neces- 
sary can be completed in reasonable time. These require- 
ments made using an adaptive-mesh grid to represent 
the simulation geometry a necessity, and the code was 
parallelized to handle the demanding computational re- 
quirements. Another desired goal was the ability to effi- 
ciently predict spectral features, such as emission lines, 
the 4000 A break and the UV continuum slope, that are 
used in observational studies of galaxies. The "polychro- 
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matic" implementation of the Monte Carlo method pre- 
sented here achieves this goal. While hydrodynamic sim- 
ulations would seem to be an ideal source for geometry- 
dependent radiative-transfer problems, this has not been 
done more than occasionally. Those cases have either con- 
cerned protostaxs_JFi^ch*r^t_al|l99J) or star-forming re- 
gions IjKurosawa et al] 120041 ). or, if they have simulated 
galaxies, have used more schematic trea tment of the hy- 
drody na mics and the radiative transfer IIBekki fc Shioval 
l2000allbl ICattaneo et all I200jtF | . To our knowledge, this 
work is the first that combines SPH hydrodynamic simu- 
lations including supernova feedback with a full radiative- 
transfer model to study the effects of dust in gala xies. Re- 
sults f rom t he simulations have b een presented in l.Tonssonl 
J2004I) and I Jonsson et all feOOd) . and additional papers 
are in preparation. 

While a number of implementations of Monte-Carlo 
radiative transfer codes are described in the literature, 
none of these are publicly available. This is in marked 
contrast to hydrodynamic codes, of which several are 
publicly available. As a service to the community, the 
author is releasing sunrise as free software. 

The organization of this paper is as follows: in sec- 
tion |5J an overview of the radiative-transfer problem 
and the Monte-Carlo method is given. Section |]3 de- 
scribes the core radiative-transfer algorithm, while sec- 
tion ^] describes the additional steps necessary to apply 
the radiative-transfer calculations to the output of hy- 
drodynamic simulations. The code is tested against pre- 
viously published results in section In section |S| im- 
plementation details are given and finally prospects for 
future improvements are given in section Q 



2 BACKGROUND 

2.1 The Radiative- Transfer Problem 

The radiative-transfer problem is the problem of calculat- 
ing the propagation of radiation through a medium which 
may emit, absorb or scatter the radiation. In the case of 
the problem of the transfer of stellar radiation through 
astrophysical dust, there are a number of simplifications 
that can be made, which greatly reduces the complexity 
of the problem. First, astronomical systems evolve slowly 
enough that the time dependence practically always can 
be ignored, and the equation to be solved is the time- 
independent radiative-transfer equation: 

k-vi v + P k v i v = p (^- + c a / Mk k')L(k')dn' 

(i) 

The dependent variable, /„, is the "specific intensity", 
defined by 



AE - I u (k,x,t)k ■ dAAQdudt. 



(2) 



In other words, it is the amount of energy per unit time 
per area per solid angle per frequency interval flowing 
in direction k. k v is the total opacity of the medium. 
The terms on the right-hand side of Equation^represent 
emission of radiation and the addition to the intensity 
from radiation scattered into the direction k from other 
directions. 



The emission from stars is fixed and does not depend 
on the local radiation field. There is emission from dust 
grains, thermal emission by the grains which are heated 
by the radiation field, but the dust grains emit mostly 
at wavelengths where stars do not and where the dust 
itself is largely optically thin. This means that the con- 
tribution from dust emission at wavelengths where stellar 
emission is important can be ignored except for a small 
wavelength interval around 5 ^tm and that the contribu- 
tion to the heating of dust grains by the emission from 
other dust grains is a second-order effect only important 
in very optically thick regions. 

Furthermore, we do not require knowledge of the full 
radiation field at all points in space; normally it is suffi- 
cient to know the radiation field that escapes to infinity, 
which is what reveals the external appearance of the ob- 
ject, and the average intensity at points inside the object, 
which is what determines the heating of the dust grains. 
It is the presence of the scattering integral in Equation 
which poses the largest difficulty. Scattering, however, is 
easily accounted for when using Monte Carlo methods. 



2.2 The Monte-Carlo Method 

The Monte-Carlo method is a way of solving equations 
by random sampling, and its usefulness for solving par- 
ticle transport problems was first recognized by Fermi, 
Ulam, and von Neumann during the days of the Manhat- 
tan Project. For an overview of the Monte-Carlo method 
as it pertains to partic le transport problems, see e.g. 
Duoree & Fralcv (2002) or, for more advanced topics, 
iLux fc Koblingerl (ll99lT) . Essentially, Monte Carlo is the 
"natural" way of solving the radiative transfer equation, 
in the sense that the photons in nature are unaware that 
they are solving a very difficult equation, they are simply 
obeying the local conditions. In the same way, solving the 
RTE using the Monte-Carlo method amounts to statis- 
tically sampling the processes that emit, scatter and ab- 
sorb photons. After tracing many such photons, statistics 
build up and the solution emerges. This local characteris- 
tic also makes the Monte-Carlo method naturally able to 
handle arbitrary geometries and complicated scattering 
characteristics 



2.3 Drawing Random Numbers 

Random numbers are at the heart of the Monte-Carlo 
method. The ability to draw random numbers with vari- 
ous probability distributions is essential. 

Given a Probability Density Function (PDF) f(x) 
for a stochastic variable X, if we generate a random num- 
ber R, uniformly distributed on [0, 1], solving the equa- 
tion 



f(x')Ax' = R 



(3) 



yields a number x with a PDF f(x). (R will be used 
throughout this paper to denote a realization of a 
stochastic variable distributed uniformly over [0, 1]. This 
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means that R in one expression is never equal to R in 
another, just that they were drawn from the same PDF.) 
As an example, let us derive the PDF of the optical 
depth at which a photon will interact. Photons propagat- 
ing through an opaque medium interact after traversing 
different optical depths — most quickly, before travers- 
ing unit optical depth, while a few make it through many. 
These optical depths are distributed as 



P(t) 



(4) 



The optical depth n at which a photon will interact with 
a medium can then be randomly generated using Equa- 
tion |3] 



- T dr = R 



n = - ln(l - R) 



InR. 



(5) 



(The last equivalence in Equation[I>]may look bizarre un- 
til it is realized that 1 — 7? has the same PDF as R\) 
Equation is the formula used to randomly draw inter- 
action optical depths in the code. 



2.4 Biasing 

It should be noted that Equation is not a unique (and, 
depending on the situation, not even necessarily the best) 
way of sampling interaction optical depths. The proba- 
bility distribution from which a quantity is drawn can 
be "biased" to sample certain parts of the distribution 
more heavily. In order to preserve a statistically correct 
result, the difference in probability must be compensated 
for by assigning a "weight" to the samples. In general, if a 
process with the probability distribution f(x) is sampled 
with another probability distribution g(x), the samples 
must be weighted by the factor w = f(x)/g(x). The only 
requirement on g(x) is that g(x) > for all x for which 
/(*) > 0. 

Biasing can be a very effective way of reducing the 
variance in the results of a Monte Carlo calculation by 
more effectively sampling the parts o f the distribution 
that are important for the end result IIDupree fc Fralevl 
I2002T ) . Ijuve la (2005) explored different biasing schemes 
in the context of radiative transfer through a spherically 
symmetric dust cloud, and found potential increases in 
efficiency by more than an order of magnitude. However, 
selecting a proper biasing requires a priori knowledge of 
the specific problem. Biasing is also the theoretical basis 
for the polychromatic algorithm described in Section l3,8l 
since it allows different wavelengths, with different mean 
free paths, to be sampled with identical rays. 



3 THE RADIATIVE-TRANSFER 
ALGORITHM 

As explained previously, the radiative-transfer problem 
will be solved through statistical sampling of the pro- 
cesses of photon emission, scattering and absorption. 

In sunrise, the dust opacity is represented on an 
adaptive grid, the characteristics of which are described 



in Section f4.3l There are many possible sources of emis- 
sion, for example collections of point sources, external dif- 
fuse radiation or emission continuously distributed on the 
same adaptive grid. In the galaxy merger simulations, the 
emission comes from the finite-sized "stellar particles" 
tracked by the SPH code. The use of an adaptive grid 
allows the representation of arbitrary geometries, limited 
only by the amount of computer time dedicated to run- 
ning the problem. (In principle, memory is also a limita- 
tion, but in practice it has been found that memory use 
is much less of a constraint.) 



3.1 Ray Tracing 

The simplest implementation of the Monte-Carlo 
radiative-transfer algorithm follows a single photon 
through the medium. This photon is emitted in a ran- 
dom direction and can then scatter and/or be absorbed. 
Eventually the photon leaves the medium in some direc- 
tion, which in general is not the direction from which the 
object is being imaged. This method is in general very 
inefficient. The efficiency can be greatly increased by cal- 
culating some probabilities analytically by having each 
ray represent a (large) number of photons, a "photon 
packet", whose intensity is proportional to the luminos- 
ity in the packet. This makes possible the use of inher- 
ently probabilistic constructs like the dust grain single- 
scattering albedo (the ratio of the scattering to the total 
cross-section) to determine the intensity of scattered ra- 
diation, rather than explicitly Monte-Carlo sampling the 
absorption and scattering processes, which is much less 
efficient. In general, analytic calculations are more effi- 
cient than explicit Monte-Carlo realization and should 
be used whenever possible. 

Furthermore, since the main object is to generate 
images of emerging radiation, it is too inefficient to sim- 
ply let rays emerge in random directions. Instead, an al- 
gorithm which efficiently estimates the radiation which 
would emerge in the directions from which the system is 
imag ed is needed. SUNRISE u ses the "Next-Event Estima- 
tor" JDupree fc Fralcv 2002), which efficiently calculates 
the flux at a point. (This m ethod is also described in 
e.g. lYusef-Zadeh et alJ 11984 .1 Using this estimator, the 
contribution to the radiation field at the observer is cal- 
culated analytically at each point of emission and scatter- 
ing. The algorithm is described in the following s ection , 
using a formalism similar to that of lGordon et alJ |2001). 
For an explanation of the symbols used, see Table 

For the ray tracing, every wavelength is treated in- 
dependently. This approach is valid since scattering by 
dust grains is an elastic process; the photon is either com- 
pletely absorbed or scattered without changing its wave- 
length. Most Monte-Carlo codes either trace rays at a set 
of discrete wavelengths, or the wavelengths of the rays 
are sampled randomly from an appropriate probability 
distribution, sunrise has until now used discrete wave- 
lengths, but here a new approach, where biasing is used 
to implement a "polychromatic" algorithm where every 
ray samples every wavelength, is presented. The descrip- 
tion that follows initially assumes that the ray tracing 
is done for one specific wavelength. The additions neces- 
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Symbol Description 



Dust opacity, or interaction cross-section per unit mass. 

Dust grain albedo (fraction that is scattered during an interaction). 

Scattering phase function asymmetry, (cos 9) . 

Luminosity of emitter e. 

Intensity (normally £ [0, 1]) of the i'th ray after the j'th interaction. 

(j = is the point of emission.) 
Total luminosity of the system. 
Total number of rays traced. 
Luminosity normalization, Ltot/JVtot- 
SPH smoothing length of the particles 
Volume of cell c. 
Path length for which the i'th ray, after the j'th interaction, passes through 

cell c. 
Luminosity associated with the i'th ray. 

Flux reaching the observer from the j'th interaction site of the i'th ray. 
Direction vector of the j'th ray after the j'th interaction. 
Direction vector towards the observer from the j 'th interaction site of the 

i'th ray. 
Distance from the the j'th interaction site of the i'th ray to the observer. 
Optical depth from the the j'th interaction site of the i'th ray to the observer. 
Optical depth from the the j'th interaction site of the i'th ray 

to the edge of the medium. 
Randomly drawn interaction optical depth of the j'th interaction of the 

i'th ray. 
Angular distribution of emitted radiation. 
Probability distribution of radiation from direction f scattering into 

direction F. 
Probability of radiation scattering an angle 9 = f • f ' (scattering phase 

function). 
Absorbed luminosity in cell c. 
Radiation intensity in cell c. 
Surface brightness in pixel p. 



9 

Ce 

Ltot 

ATtot 

n 

h 

V c 



h ■ 
j£,obs 



_obs 



-MM') 

<& s (cOS0) 

A c 

Jc 



Table 1. Symbols used in the description of the radiative-transfer algorithm and their meaning. There is an implicit dependence 
on wavelength in all these quantities. 



sary for a polychromatic algorithm are then described in 
Section [ 



3.2 Dust Properties 

In order to perform the ray tracing, the properties of 
the opaque medium must be specified. The necessary 
quantities are k, the mass opacity of the dust grains; 
a, the single-scattering albedo; and finally the scatter- 
ing phase function, the angular distribution of the scat- 
tered photons. While sunrise is capable of using an ar- 
bitrary phase function, the one currently used is the pop- 
ular ^albejtoi^uestJOTiabl^jcCTn^c^jas poi nted o ut by 
e.g. lDrainei200a)lHenvev fc Greensteinl (HG. Il94ll) func- 
tional form, 



$ s (cos 8) — 



9' 



(6) 



4tt(1 + 5 2 -2 5 cos(9) 3 / 2 ' 

where 8 is the scattering angle and g — (cos 8) , the phase 
function asymmetry, parameterizes the degree to which 
the scattering is isotropic or mostly forward/backward. 
The HG function has the advantage that it can be ana- 
lytically inverted. 

In order to calculate the absorption and scattering 
of light by the dust, all that is needed are the three quan- 
tities k, a, and g, as a function of wavelength. Knowledge 



about the detailed grain size distribution and composi- 
tion giving rise to these quantities is not necessary. How- 
ever, if the spectrum of infrared radiation emitted by the 
dust is to be calculated, the situation is different. When 
this capability, which is a planned upgrade, is added to 
SUNRISE, these details will be necessary. 

3.3 Emission 

The luminosity associated with ray i at emission is 

Li = Iifln , (7) 

where Iij is the intensity of the ray, identifying its statis- 
tical weight, and n = L to t/Ntot is an overall normaliza- 
tion (common for all rays). When rays are emitted, the 
probability of emission at a certain point and in a cer- 
tain direction is normally equal to the actual probability 
of emission of photons, and their initial intensity is unity. 
(As explained in Section [2. 41 this is merely a choice.) 

The reason for the separate normalizing factor n is 
that it is often desirable to be able to delay the final 
normalization to until after all rays have been traced, 
while It must be known at the time the ray is created. 
This way, the total number of rays A^ot does not have to 
be known in advance, which can be the case for example if 
the ray tracing is being done on several CPUs in parallel. 
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In the case of several sources of emission, like a col- 
lection of particles or many grid cells, the emitter e from 
which the ray is emitted is drawn randomly by finding e 
such that 



J2 P e' < R <Y1 Pe ' 



(8) 



where P e — L e /L to t is the probability of emission from 
emitter e. 

Once the source of emission has been determined, 
the originating position x e is determined. If the emitter 
is a finite-sized particle from the SPH simulations, the ra- 
dius of the point of emission is determined based on the 
density profile represented by the SPH smoothing ker- 
nel. In order to avoid resorting to numerically solving for 
the radius, a polynomial approximation of the probabil- 
ity distribution resulting from the SPH smoothing kernel 
is used. The simplest polynomial representation that has 
a mass which goes to at small radii and a density which 
goes to at twice the smoothing length h is 



P(r < ah) 



2, 

-a [a 



(a - 3)/4 , < a < 2 



0) 



which corresponds to the density distribution 

rt«) = 4(?-i)- (io) 



The radius of emission is then determined by solving the 
cubic equation P(r/h) = R. 

If the source of emission is a grid cell, the position 
within the cell is drawn from a uniform random distribu- 
tion: 



— •^min,c T [X 



max,c 



^min,c 



(Rx + Ry + Rz) (11) 



where i m i„, c and x ma x,c are the lower and upper bound- 
aries of the cell. (Remember that the three instances of 
R denote three different random numbers.) 

With the point of emission determined, the direct 
flux that would result from the emission, if the ray was 
emitted in the direction of the observer, is calculated: 



Fi, = nlifie 



D ° a $e(fc° b o s : 



«*?,„' 



(12) 



obs 



Tj"o" is the optical depth between the point of emission 
and the observer, "J> e is the angular distribution of emit- 
ted radiation and k° q s the direction vector from the site 
of emission to the observer. (In the case of isotropic emis- 
sion, $ = l/(47r).) di t o is the distance from the site of 
emission to the observer. This calculation is repeated to 
calculate the contribution in all "cameras" . 

Finally, a specific direction of propagation fe;,o for 
the ray is randomly drawn from the angular distribution 
of the emitted radiation, <3? e , defined as 



!>,: 



! 



d/ 
dfi 



such that 



<J> e dfi = 



AI_ 



dtt- 



(13) 
(14) 



Using Equation^ it is possible to draw random directions 
from this distribution. In the case of isotropic emission, 



the procedure is 

9 = arccos(2i? - 1) , 
(j) = 2-kR , 
ki,o — (sin 9 cos <j>, sin 9 sin <f>, cos ( 



(15) 



The code also includes other types of emitters, such as 
point sources, and other angular distributions of the emit- 
ted rays, like collimated beams. Arbitrary angular distri- 
butions are easily added. 



3.4 Ray Propagation 

The ray is now propagated in the direction fe;,o- The 
propagation is done from one cell to another, keeping 
track of the optical depth traversed by the ray. At each 
step, the optical depth is increased by At = p c nAl, 
where p c is the density of dust in the cell, k is the opacity 
of the dust, and A£ is the path length traversed by the 
ray inside the cell. 

If the medium traversed is optically thin, most rays 
would leave the simulation medium without interacting 
and not contribute to the scattered flux. To increase the 
efficiency of the calculation of the scattered flux, SUN- 
RISE, like most other Monte -Carlo codes, uses the con - 
cept of "forced scattering" JCashwell fc Everett! Il959f) . 
in which every ray is forced to contribute to the scat- 
tered flux. In the "forced scattering" scenario, the to- 
tal optical depth rf from the point of emission x e to 
the edge of the medium in the direction of propaga- 
tion is first calculated (by tracing the ray to the edge 
of the grid) . The ray is then split up into two parts. One 
part, /i,oe _Tio , will leave the medium without interac- 
tion, while Ji,o ( 1 — e~ T '-° ) will interact somewhere along 
the path. The optical depth of this interaction, which is 
in the range [0, t^ ], is drawn randomly using the formula 



Tift 



■ In 1 - R 1 



(16) 



which is a variant of Equation [2] obtained by restricting 
the range of optical depths to [0, rf ] an d renormaliz- 
ing the distribution. The part of the ray that leaves the 
medium is dropped, as the flux resulting from direct ra- 
diation already has been taken into account with Equa- 
tion^] The part of the ray that does interact will have 
an intensity after the interaction of 



h,l — lifiO, 1 1 



(17) 



where a is the dust grain albedo. The luminosity absorbed 
in the grid cell where the interaction takes place is 



At,i = nIi,o(l - a) (l 



(18) 



The part of the ray left after the interaction is scattered 
into a new direction by the dust grain. Analogously to 
Equation 1121 the flux resulting from the part of the ray 
which would be scattered towards the observer and which 
would not interact on its way there, is 

F iA = n/i,ie-^*.(fci,o, fc?,i s ) J-, (19) 

where <J? s (fe,fe') is the scattering phase function, i.e. the 
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angular distribution dl/dQ for scattering of rays from 
direction k into direction k . In most cases, the phase 
function will depend only on the angle between the two 
directions (exceptions to this would be e.g. non-spherical 
dust grains aligned by a magnetic field), in which case 
$ s (fe, k') = <!> s (fe ■ k'). As mentioned in Section l3~2l SUN- 
RISE currently uses the Henyey-Greenstein function, but 
is capable of using arbitrary ph ase functions. (Note that 
the equation equivalent to lliJI in lGordon et alJl200ll their 
equation 6, has an extraneous Air [K. Gordon 2004, pri- 
vate communication].) 

After calculating the intensity that would have made 
it to the observer, had the ray been scattered in that 
direction, a random scattering angle is drawn from the 
scattering phase functio n. In the ca se of the HG phase 
function, the formula is <IWitdll977ab 



cost/ = 



1 
2%R 



i i ^ 
1+3 - 



g + 2gR 



(20) 

(21) 



(The arbitrary reference axis for the azimuthal angle <j> is 
taken to be the z axis.) The ray is then rotated by this 
angle, and this becomes the new direction of propagation. 
Depending on the problem at hand, the number of 
scatterings that are forced may range from up to any 
number. Normally, only the first scattering is forced, but 
if one is interested in the effects of higher-order scat- 
tering, a larger number may be motivated. If the j'th 
scattering is not forced, an interaction optical depth is 
drawn according to equation |S] The ray is then propa- 
gated through the medium until it either leaves or reaches 
the interaction optical depth. If it reaches the interaction 
optical depth, it is scattered. The intensity of the ray 
after the interaction is then 



lij = alij-i , 



while the absorbed luminosity is 

Aij — (Ii,j — Iij-i)n = (1 — a)nli. 



i-i 



(22) 



(23) 



The flux resulting from the intensity scattered in the di- 
rection of the observer is 



fit — nJ. 



"Mfc; 






d 2 ■ 
a '.j 



(24) 



This process is repeated until the ray either leaves 
the volume or until the intensity of the ray drops below 
some threshold, 7 m in, set by the user. To avoid expend- 
ing computational resources tracking rays with very low 
intensity that will not contribute significantly to the re- 
sults, the ray is dropped once its intensity drops below 
Imin- However, to avoid violating energy conservation, 
this cannot be done in all cases. Instead, the ray is given 
some probability Prr of survival, and a random number 
is drawn to see if the ray survives. If it does, its intensity 
is increased by a factor 1/Prr and the ray continues to 
be tracked. If not, the ray is terminated. This scheme, 
known as "Russian roulette", ensures that energy con- 
servation, on average, is preserved. 



3.5 Multiple Scattering Components 

It is possible to define an arbitrary number of scattering 
components in every grid cell, for example if there are two 
different types of dust with radically different scattering 
properties, distributed differently. In this case, the optical 
depths in the formulae above refer to the sum of the op- 
tical depths of all the components. When an interaction 
takes place, the component responsible for the scattering 
event is drawn randomly with a probability proportional 
to the opacities of the components in the grid cell where 
the interaction occurs. The scattering is then performed 
using the albedo and scattering phase function of this 
component. 

It should again be pointed out that this procedure 
applies to the transfer of radiation through the medium, 
for which only aggregate quantities is necessary, and not 
to the determination of grain temperatures. For exam- 
ple, if the medium contains both dust with Milky- Way- 
like and Small-Magellanic-Cloud-like properties, both of 
which represent a distribution of grain sizes and composi- 
tions, this procedure is used when calculating the transfer 
of radiation through this medium. Only if the grain tem- 
perature distributions, which are dependent on the full 
set of grain sizes and compositions, is to be calculated is 
it necessary to consider each grain size separately. 



3.6 Output Images 

The rays that are calculated to make it to the observer 
are projected through a virtual "pinhole camera" onto an 
image plane. These cameras can be placed at arbitrary 
points. (However, when cameras are placed at a position 
where emission or scattering can occur, the noise in the 
images will increase drastically due to the infinitely large 
contributions from events occurring infinitesimally close 
to the camera position. This p roblem is known as th e 
"infinite variance catastrophe" JDupree fc Fra lev 2002), 
but is not likely to occur in astronomical simulations.) 
Each camera has a specified field of view and image array 
size. 

The surface brightness of pixel p in the image is then 
calculated as 



Zjp — 



i 'nix 



(25) 



where the sum over i (ray number) and j (interaction 
number) only includes those Fi.j whose point of origin 
project onto pixel p, and O pix is the surface angle sub- 
tended by the pixel for the observer. If the projected ray 
is outside the field of view of the camera, the light is lost. 



3.7 Absorbed Luminosity 

The absorbed luminosity in grid cell c is calculated as 

Ac = Y,Ai,j, (26) 

i,i 

where the sum over i and j only includes those ray inter- 
actions that occur in cell c. The total absorbed luminosity 
in a grid cell equals the total luminosity reradiated by the 
dust, by energy conservation. However, if the dust grain 
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temperature distribution and the SED of the dust emis- 
sion is to be calculated self-consistently, the (wavelength- 
dependent) radiation field in the cell must be determined 
( Guh athakurta fc Drainell989l ) . If the absorbed luminos- 
ity is known, the radiation field J c in the cell can be cal- 
culated as 



MX) = 



A c 



47Tp c « abs V c ' 



(27) 



where K a b s is the absorption opacity of the dust and V c 
is the volume of the cell. Because only absorption events 
contribute to the signal, this method suffers from large 
Monte-Carlo noise in regions where the number of inter- 
actions are few, for example in highly refined cells with 
small volume. In fact, because the radiative-transfer algo- 
rithm used by sunrise is so much more efficient at getting 
signal to the cameras than the simplest Monte-Carlo im- 
plementation, fewer rays need to be traced. This means 
that, since each ray interacts with the medium at most a 
few times, the number of absorption events determining 
A(X) in eauation l27l is small and the qu antity noisy. SUN- 
rise uses another met hod, described bv lLucvl II1999T) and 
iNiccolini et al.l a2003l) . that takes advantage of the fact 
that, physically, the radiation intensity is determined by 
the number of rays (photons) traversing a volume, re- 
gardless of the probability of absorption. In this scheme, 



J2i d &li,j,c.nli, 



(28) 



where A^; JiC is the path length during which the i'th 
ray, after the j'th interaction, passes through cell c. Since 
many more rays pass through a given cell than are ab- 
sorbed in it, this method has superior accuracy. The only 
complication is in the case of forced scattering. In this 
case the ray intensity in the cells traversed before the 
forced scattering takes place is J^o, but the part of the 
ray that leaves the medium without interaction also has 
to be taken into account. That part of the ray has lower 
intensity; it gives a contribution to J c corresponding to 
Ii,o&~ Ti '° in the cells traversed after the interaction point. 



3.8 Polychromatic Ray Tracing 

Looking back at the preceding sections, it is possible to 
identify the points where wavelength dependence poses 
conceptual problems to a procedure where all wave- 
lengths are included in every ray. Clearly, emission of 
rays and calculation of the flux reaching the observer ei- 
ther directly (Equation I12H or from a point of scattering 
( Equation 1191 pose no problems: these formulae are an- 
alytic calculations that can be performed for any num- 
ber of wavelengths simply by replacing the quantities h, 
t°j b and a with vectors of numbers. Polychromatic cal- 
culations of t he direct flux ha s already been done in the 
SKIRT code feaes et al.l2005i) . The problems arise where 
interaction optical depths and scattering directions are 
sampled from the appropriate probability distributions, 
because these probability distributions depend on wave- 
length. For example, rays of shorter wavelength will tend 
to travel shorter distance before interacting, since the 



dust opacity generally increases towards shorter wave- 
lengths. This means an interaction point can only be 
drawn in a statistically correct way for one wavelength 
at a time, and the same objection applies to the scat- 
tering angle drawn using eouation l2UI However, the abil- 
ity to use biased distributions opens up the possibility to 
compensate for the fact that the probability distributions 
will only be correct for one wavelength. (This is known 
as "path stretching".) The proper way of doing this will 
now be examined. 

As was derived in Section f2.3l the probability distri- 
bution function of where a ray interacts with the medium 



dP(r) 



~ T dr. 



(29) 



Suppose an interaction optical depth r rG f is drawn for 
some reference wavelength A rc f. The probability of an- 
other wavelength A interacting at the same point is then 



dP[r(A)] = e- T(A) dr(A) 



r(A) 

T re { 



dr r , 



(30) 

The necessary biasing factor wx is the ratio of the prob- 
abilities at wavelengths A and A re f: 



w x 



PJr(X)] 

P [Trrf] 



8 f-r(A) 



r(A) 



T rc f 



(31) 



To compensate for the biased probability distribution, 
the intensity of the ray at different wavelengths at the 
point of interaction should be multiplied by the weight- 
ing factor wx before calculating scattered or absorbed 
luminosity. 

In cases where forced scattering is used, the probabil- 
ity distribution from which interaction points are drawn 
is different, and so is also the weighting factor. The cor- 
rect Wx when forced scattering is used it is 



W\ — e 



T rof -T(A) 



~(X) 



T re f 



1 



1 



"(>■) 



(32) 



Finally, the biased distribution of scattering angles 
must be accounted for. Compared to the optical depths, 
this is quite straightforward: the probability of scattering 
into a certain direction is given by the scattering phase 
function & s (9), so if a scattering angle 6 is drawn at the 
reference wavelength the weighting factor which should 
be applied to the ray intensity after scattering will be 



Wx 



M8,X) 

<MMrcf) 



(33) 



Energy conservation, in a statistical sense, must be 
maintained in the polychromatic calculation; energy flux 
is the product of probability and ray intensity, and any 
biasing scheme simply trades probabilities for intensities. 

The possibility of calculating all wavelengths simul- 

I 1 J" 

taneously was noted bv lJuvelal (2005). who argued that it 

would not be advantageous since the opacity is a strong 
function of wavelength and the large bias factors nec- 
essary probably would result in increased errors. It is 
true that the errors for a fixed number of rays probably 
would increase for wavelengths where the dust opacity is 
very different from what it is at the reference wavelength, 
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but the differential errors between different wavelengths 
are minimized. This is clearly illustrated in the example 
spectra calculated in Section 15.41 Because every wave- 
length is uncorrelated in the monochromatic calculation, 
the spectral shape becomes very noisy in regions of low 
signal-to-noise, but this is not the case with the polychro- 
matic calculation. (The use of "correlated Monte Carlo" 
for perturbation analysis builds on the same principle, 
the stochastic effects are minimized by correlating the 
random walks in the perturbed and unperturbed cases.) 
The increased noise at wavelengths far away from the 
reference wavelength is alleviated by the fact that with 
the polychromatic algorithm more rays can contribute to 
each wavelength. This is because the computational costs 
of tracing rays in sunrise is dominated by propagating 
rays from cell to cell during the random walk. As long as 
the vector operations for doing the calculation at many 
wavelengths is not the dominant computational cost, the 
extra wavelengths are obtained at little cost. 

A great advantage of the polychromatic algorithm 
over the approach using a set of discrete wavelengths 
is that the radiative-transfer problem is truly solved for 
every wavelength; spectral features present in the stel- 
lar emission are predicted properly, and the differential 
attenuation between lines and continuum is accurately 
treated. (This can also be seen in the example spectra cal- 
culated in Section 15. 41 ) This is important for predicting 
images and spectra of galaxies, the particular application 
for which sunrise was developed, where differential ex- 
tinction between different stellar populations can signifi- 
cantly affect spectral features like the Balmer absorption 
lines. 

Another advantage is that it facilitates the inclu- 
sion of kinematic effects into the radiative transfer. Kine- 
matic effects are already included in the SKIRT code 

I II 3\ 

(Bacs ct al. 2003), but then only a small range in wave- 
lengths, over which the dust characteristics is assumed 
to be constant, is included in the calculation. Using 
the polychromatic algorithm, the entire spectrum can be 
propagated as the accumulated Doppler shift is tracked. 
When the ray is finally projected on to the camera, this 
Doppler shift can be included yielding a full spectrum 
including Doppler broadening. 



3.9 Uncertainties 

Because of the stochastic nature of the Monte-Carlo 
method, results are subject to random sampling error. If 
this error can be evaluated at runtime, the number of rays 
used can be adapted to get the error below some toler- 
ance iJGordon et al.l20Qj) . and even if this is not possible, 
knowing the uncertainty in the results is important. 

Estimating the uncertainty is conceptually straight- 
forward. The quantity of interest (the flux in a pixel, for 
example) is the sum of Nv samples of a random variable 
V. The variance in the estimate of the sum of Nv sam- 
ples of a random variable is N v times the variance of V, 
which in turn is estimated as 



o 10 




2 



1 V^ 2 1 v^ 



n 



(34) 



10 5 10 5 10 7 
Number of rays 



Figure 1. The fractional errors (standard deviation/ value) 
of the pixel values in the image shown in Figure[I]for different 
numbers of rays, showing the convergence of the results. The 
uncertainties were estimated both using equation 1341 (red di- 
amonds) and as the empirically determined variance from 30 
runs (blue triangles). The error bars show the la logarithmic 
spread in the distribution of pixels, and the solid line indi- 
cates the theoretically expected y/N convergence. For small 
numbers of rays, equation l34l severelv undcrprcdicts the actual 
variance. (The standard deviation estimated with equation l34l 
is bounded by the quantity itself in cases where only one ray 
contributed to the result.) For increasing numbers of rays, the 
estimates converge towards the empirically determined values, 
but does so systematically from below. There is a tendency 
for the results obtained with the "next-event estimator" to 
be dominated by rare, large contributions, implying that un- 
less the problem phase space is well sampled, the estimated 
uncertainty will always be too small. 



However, there is an important subtlety here. A "sam- 
pling" of the flux in a pixel consists of the shooting of 
one ray. Vi is the total contribution made by the ray, so 
if one ray contributes flux in the pixel both directly and 
through subsequent scatterings, the sum of all these con- 
tributions is Vi. This is important, because Vi is squared. 
Squaring the contributions from direct, single-scattered, 
etc., light separately will lead to a systematic underesti- 
mate of the variance. It should be pointed out that, for 
large Nv , the estimated variance is insensitive to whether 
Nv is taken to be the total number of rays or just the 
number of rays with nonzero contributions to the pixel. 
(In general, most rays will not contribute to the flux in a 
given pixel.) 

This formula for estimating the uncertainty has been 
tested on the pixel-by-pixel brightnes s in one of the 
cases o f the clumpy scattering medium of Wit t fc Gordon! 
( 1996) against which the code is tested in Section 15.31 
This test is shown in Figure For different numbers 
of rays, the estimated variances in the pixel values from 
equation 1341 were compared against the empirically de- 
termined variances from 30 runs with different random 
numbers. For small numbers of rays, there are many pix- 
els which have one or no contribution, and the variance 
is severely underestimated. As the number of rays in- 
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creases, the variances estimated using equation 1341 con- 
verge toward the empirically determined values, but does 
so from below. This emphasizes an important fact about 
variance estimates using the next-event estimator: the re- 
sults tend to be dominated by rare, large contributions, 
so unless the problem phase space has been thoroughly 
sampled, the true variance will be underestimated. (In 
fact, if scatterings can take place close to the camera, the 
true variance is not even finite, as previously mentioned.) 
As seen above, it is necessary to save not only the 
quantity itself but also the sum of the squared contri- 
butions to determine the uncertainty. This increases the 
amount of data to keep track of by a factor of two, and 
makes the parallelization of the algorithm more compli- 
cated. Because of this, the uncertainty is normally not 
estimated in the galaxy merger simulations (The outputs 
from one run are already frequently larger than 1 Gb) . 



written, additional steps are necessary. These extra steps 
are described in this section. 

The general procedure is as follows: first, hydrody- 
namic simulations are used to generate the geometry of 
the problem, e.g. merging galaxies. A number of snap- 
shots at different time steps are saved, and for each of 
these a series of preparatory steps is performed. First, 
a stellar population synthesis model is used to calcu- 
late the SEDs of the emitting sources. Second, the adap- 
tive grid needed for the radiative-transfer calculations is 
generated. Third, the radiative-transfer calculations are 
done. (Work is currently being done to integrate the poly- 
chromatic version of SUNRISE into this framework. The 
hydrodynamic simulations processed so far have used to 
the monochromatic version of sunrise.) Finally, a post- 
processing step is done, where the full SED is calculated 
by interpolating over wavelength. 



3.10 Clumpy Dust Distributions 

Many authors have studied the effects of clumpy distribu- 
tions of dust JWitt fc Gordonlll996t f2000: Bia nchi et all 
2000) and reached the conclusion that it can profoundly 
affect the attenuation of radiation for a given mass of 
dust. The general tendency is for clumping to decrease at- 
tenuation and reddening by allowing radiation to escape 
through optically-thin lines of sight, unless the sources 
themselves are embedded in the clumps, in which case 
the attenuation can be higher than the corresponding 
homogenous scenario. 

In sunrise, dust is assumed to be uniformly dis- 
tributed within each grid cell. There is no additional 
clumping assumed beyond what is actually present in 
the adaptive grid. In our hydrodynamic galaxy merger 
simulations there are large-scale inhomogeneities, such 
as spiral structure, present, but the resolution (and the 
physics incorporated, for that matter) is too coarse to 
resolve e.g. individual molecular clouds. The adaptive 
grid structure could be used to put in artificial clump- 
ing on yet smaller scales, but the computational require- 
ments would be prohibitive. One solution would be to 
incorporate clumping through a sub-grid analytical ap- 
proximation, such as the "mega-grains approximation" 
<lNeufeldlll99l iHobson fc PadmarJlilffll IVarosi fc Dwekl 
|l999j); large clumps can be treated as enormous dust 
grains, with an effective albedo and scattering function. 
Clumps of dust can then be added within the existing 
framework as just another source of scattering. If the 
sources of emission are located within the clumps, how- 
ever, this approach will not work. This can instead be ac- 
commodated by changing the characteristics of the emit- 
ted radiation before injecting the rays into the grid. 



4 AUXILIARY CALCULATIONS - 
RADIATIVE TRANSFER IN 
HYDRODYNAMIC SIMULATIONS 

In section [!|1 the core radiative-transfer algorithm of SUN- 
RISE was described. However, in order to generate images 
and spectral energy distributions from hydrodynamic N- 
body simulations, the purpose for which the code was 



4.1 The Hydrodynamic Simulations 

The hydrodynamic s imulations have been described 
in detail elsew here lICoxl |2004 Ijonsson et alJ l200rj : 
ICox et all2005l) . In particular, the specific galaxy models 
which are bein g used for the galaxy merger simulations 
are described in Jon sson et al.l l|2006|) and will not be dis- 
cussed here, but in order to provide a context and define 
the quantities used in the radiative-transfer calculations, 
a brief overview of the technique used will be given. 

The si mulations are done u sing GADGET 
(ISpringel et alJ l200lt ISpringcl 2005), a Lagrangian 
Smooth Particle Hydrodynamics (SPH) code. The 
galaxies are initially modeled as a disk of stars and 
gas, a stellar bulge, and a dark-matter halo. The stars 
and dark-matter particles are collisionless and only 
feel the force of gravity. The gas particles are also 
subject to hydrodynamic forces. A collisionless particle 
i is characterized by its mass rrn and its gravitational 
softening length n. A gas particle has, in addition, an 
associated SPH smoothing length, hi, which indicates 
the size of the region over which the hydrodynamic 
quantities associated with the particle are averaged. 
The smoothing length is determined by the distance to 
neighboring gas particles, i.e. by the gas density, such 
that the resolution is higher in high-density regions and 
lower in low-density regions where the particles represent 
a very smoothed-out gas density field. 

During the simulation, the star-formation rate of 
each gas particle is calculated according to a "Schmidt 
law", 



dp, 
di 



1.5 
^ rgas 



(35) 



As stars form, gas is transformed into collisionless mat- 
ter. In the simulations, this is implemented in a stochas- 
tic sense (Sprinecl & Hcrnauist 2003) in which each gas 
particle spawns a number of stellar particles with a prob- 
ability consistent with the calculated star-formation rate. 
These "new star" particles represent associations of single 
stellar populations, though their mass, typically 10 6 M©, 
is larger than most observed young star clusters. In ad- 
dition to the quantities associated with all collisionless 
particles mentioned earlier, new star particles are also 
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characterized by a formation time £f,i and a metallicity 
Zi, which is the metallicity of the gas particle from which 
it is spawned. 

Associated with star formation is supernova feed- 
back, whereby energy from supernovae is deposited into 
the interstellar medium. This energy heats and pressur- 
izes the gas and stabilizes it from further gravitational 
collapse. Including feedback is crucial for the stability of 
the simulations, but it is a complicated subject and many 
different approaches to implementing it exist. Here, it will 
only be mentioned that our simulations contain feedback 
and that a significant effort has gone into c onstraining 
this part of the simulations IICox et al] 120051) . Feedback 
from AGN could also affect the g as. This has been in - 
cluded in GADGET simulations JSpringel et alJl2005h . 
and these simulations can also be processed by sunrise. 

A third consequence of star formation is chemical 
enrichment. Since the goal is to simulate the effects of 
dust, tracking metal production is naturally a topic of 
interest as the amount of metals available will affect the 
amount of dust present. In the simulations, metal produc- 
tion is implemented using an "instantaneous-recycling" 
scheme: stars formed are assumed to instantly become 
supernovae, and the metals produced are put back into 
the gas phase of the particle. Every gas particle is, in ad- 
dition to the quantities mentioned earlier, characterized 
by a metallicity Zi. This scheme, while simple, has sev- 
eral drawbacks: First, metals do not diffuse from the gas 
particle from which they were made, and if this particle is 
completely turned into stars, all the metals are incorpo- 
rated into the stars and lost from the ISM. Second, while 
metals are recycled in supernovae, gas is not. In reality, 
a stellar population returns a non-negligible fraction of 
its mass to the gas phase, due to supernovae and stellar 
winds. Third, while instantaneous recycling may be a rea- 
sonable approximation for Type II supernovae, it is surely 
not for Type la supernovae which are believed to explode 
at least several hundred million years after the stars are 
formed. Improved models f or metal enrichment in G AD- 
GET have been developed llScannapieco et aDteOOST) . but 
our simulations have so far not used these. 



4.2 Calculation of Stellar SEDs 

After the hydro simulations have been completed, the 
first step is to calculate an SED for the stellar particles 
in each simulation sna pshot. In this work, t he SEDs used 
are from Starburst99 ilLeitherer et al]ll999l) but are sub- 
sampled to 510 wavelength points in order to minimize 
file size. (Note that this is not the much smaller number 
of wavelengths for which the monochromatic radiative- 
transfer calculations are done.) The calculation of the 
stellar SED is trivial since the assumption is that stellar 
particles represent single stellar populations, so one sim- 
ply has to pick an SED with the age and metallicity of 
the particle. 

For stars present at the start of the simulation, i.e. 
the disk and bulge components of the merging galax- 
ies, the star-formation history and metallicity distribu- 
tion must be specified as input parameters. Typically, 
the bulge stars are assumed to have formed in an instan- 



taneous burst a relatively long time ago, while the disk 
has had an exponentially declining star-formation rate 
starting at the time of bulge formation and leading up to 
the start of the simulation, but any choice can be made. 
Formation times consistent with these assumptions are 
then drawn randomly for the individual particles. 

The physical size of the region over which the par- 
ticle luminosity is distributed is set to some fixed value. 
Typically, the gravitational softening length n is used. 



4.3 Creating the Adaptive Grid 

The amount of dust is based on the amount of met- 
als present in the galaxy simulations. As was de- 
scribed in section |3] the ray tracing is done on a grid, 
so it is necessary to transform the density field de- 
scribed by the particles onto a grid. In order to be 
able to resolve the small high-density regions in the 
simulations, while still covering the large volume over 
which the interaction takes place, this gr i d is adaptive 
IWolf et allf999l:lKurosawa fc HillierlgOOll ; [Harries et all 
12004 IStamatellos fc Whitwortbll2005l) . The grid struc- 
ture is based on a regular Cartesian grid, in which grid 
cells can be recursively refined by subdividing them into 
2 3 subcells, leading to an octtree-like memory structure. 

Grid construction proceeds as follows. First, the base 
grid is created. This grid covers the entire extent of the 
geometry to be simulated, and typically has 1 3 grid cells. 
Each of these cells is then recursively subdivided until 
l c < iniiii(hi)c, i.e. the cell is smaller than the sizes of all 
particles contained within it. This strategy uses the infor- 
mation contained in the SPH smoothing length to deter- 
mine the smallest scale which possibly contains structure. 
The constant c is a fudge factor, typically chosen to be 
2, since this is roughly the region over which the particle 
contributes density. (The SPH smoothing kernel extends 
to 2hi.) 

Once the refinement step is completed, the mass of 
metals contained in the particles is projected onto the 
grid cells using the SPH smoothing kernel as the radial 
density profil e. For the projection, the piecewise polyno- 
mial kernel in lHernauist fc Kata (119891) is used. Perform- 
ing this three-dimensional integral is time-consuming, 
and therefore a tabulated version is used. In the cases 
where the particle is much larger than the grid cell, the 
density associated with the particle is assumed to be con- 
stant over the extent of the grid cell, eliminating the need 
for the integration. 

At this point, the grid describes the spatial distribu- 
tion in the simulations but not necessarily in the most 
efficient way. Unlike in, for example, an adaptive-mesh 
hydrodynamics code, the necessary size of the grid cells 
does not depend on the magnitude of the density, only on 
its inhomogeneity. There is no size scale that has to be 
resolved, as long as the problem is described accurately. 
(This is not true for the determination of the radiation 
intensity, since this can obviously change even if the dust 
density is perfectly uniform. Unfortunately, there is no 
local criterion for determining the resolution required to 
resolve the radiation field without actually solving the 
radiative-transfer problem.) Given this fact, it is desir- 
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able that cells whose quantities are sufficiently uniform 
be unified into one larger cell, as the ray tracing takes 
approximately constant time per grid cell traversed. Sub- 
cells are unified and replaced with their "parent" cell as 
long as the following criterion is fulfilled: 
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where the average ((•••)) and standard deviations (a) are 
calculated over the subcells, V^ r id is the volume of the 
entire grid, N is the number of Monte-Carlo rays to be 
traced, and ft is an opacity (dr/dx). Two different kinds 
of criteria are used: The first is a relative one; the stan- 
dard deviation of the quantity over the subcells divided 
by the average quantity (the value the unified cell would 
have) must be less than a specified tolerance (tol mc t) 
for unification to be allowed. This ensures that inhomo- 
geneities are resolved. The second criterion is an abso- 
lute one; the standard deviation of the subcells must be 
smaller than some value. The idea is that if the difference 
in the quantity resulting from unification is so small that 
not even one Monte-Carlo ray will be affected by such 
a change, we can unify the cells regardless of how large 
the fractional deviation is. This ensures resources are not 
wasted on resolving regions that are so sparse they will 
not matter to the results anyway. 

If we have N rays traversing the volume V^ r id, 
the number traversing a subvolume v will be n = 
N(v/V gT id) 2 , if we assume a uniform and isotropic dis- 
tribution of rays. The number of rays that will interact 
in the subvolume is 
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interestingly independent of v. From this one can con- 
clude that a deviation in mass can be accepted without 
affecting a single scattering if 
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grid 



(38) 



Since this result is valid for a uniform and isotropic dis- 
tribution of rays, which is a questionable assumption, it 
is obviously not to be trusted in detail. It does however 
provide a scaling, and to ensure cells are not excessively 
unified, N is usually taken to be at least 10 x the number 
of Monte-Carlo rays actually used. 

The grids used for the galaxy merger simulations 
are constructed with tol mot = 0.1, N — 10 7 and 
k = 3x IO -5 kpc 2 / M©. With these parameters, the grids 
typically have 30k to 100k cells and cover a volume of 
(200 kpc) 3 , with a maximum of 10 subdivisions. This is 
the equivalent resolution of a 10240 3 uniform grid. Given 
the run-time scaling in Section 16.21 the use of an adap- 
tive grid allows the simulations to complete in 2% of the 
time that would be required for a uniform grid. Clearly 
this dynamic range, which matches the dynamic range in 
the hydro simulations themselves, would not be possible 
without an adaptive grid. 

It is also worth mentioning that although this de- 
scription was based on building the grid from the in- 
formation in the hydro simulations, the code includes a 
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Table 2. The wavelengths for which the radiative-transfer 
calculation is done. 



generic "grid factory" class which can be used to build 
the grid using density fields from any source. In addition, 
the grid can contain diffuse emission and the refinement 
also take the inhomogeneities in emissivity into account. 



4.4 The Radiative- Transfer Stage 

The radiative-transfer algorithm as described in detail in 
sectionals repeated for a number of wavelengths in order 
to get an idea of the appearance of the system at all ob- 
served wavelengths, and to determine the total amount of 
energy absorbed by the dust. The assumption is then that 
since the dust properties change smoothly with wave- 
length, the results can be interpolated over wavelength 
to give a full SED. 

The radiative-transfer stage is done in two phases. 
First, a run without considering any dust effects is done. 
The entire SED (all wavelengths at once) is simply propa- 
gated from randomly drawn emission sites directly to the 
cameras using Equation 1121 giving images of the object. 
This is possible since without dust there is no wavelength- 
dependent part to this process. These "dust-free" images 
serve as a basis for the later interpolation over wave- 
length. Second, the full radiative transfer including scat- 
tering and absorption is done for a number of wavelengths 
from the far-UV to the near-IR. This range includes prac- 
tically all the energy emitted by stars, so energy outside 
of this range is ignored. The number of wavelengths used 
is a trade-off between accuracy and computation time. 
For the galaxy simulations, 19 wavelengths, shown in Ta- 
ble |21 have been used. They were determined through a 
fitting procedure where the error in the fraction of the lu- 
minosity absorbed by dust was minimized against a test 
run using 100 wavelengths, and are particularly dense 
around the 2200 A feature in the Milky- Way dust extinc- 
tion curve to ensure that this region of the spectrum is 
calculated accurately. 

The assumption of smoothness as a function of wave- 
length is not valid if one considers line emission such as 
Ha. While the dust properties do not change over the 
line wavelength, the emission properties change abruptly; 
emission lines come from star-forming regions, which are 
generally more deeply embedded within the dust than 
are the sources contributing to the continuum around 
the line. This means that emission lines can have an at- 
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tenuation many times larger than the continuum around 
them, and as a consequence their dust dependence must 
be calculated by specifically doing the radiative transfer 
in the line. Currently, the hydrogen lines Ha and H/9 are 
included, as these are frequently used as star-formation 
and dust indicators. 



the hydrodynamic simulations. The quantities imaged are 
mass density of stars, gas and metals, star-formation rate, 
bolometric luminosity, mass- and luminosity-weighted 
stellar age, and gas temperature. These images make it 
possible to correlate the emerging radiation and the dust 
effects with physical quantities of the system on a pixel- 
by-pixel basis. 



4.5 Post-processing 

After the radiative-transfer calculation has been done, 
the results are interpolated to give an SED with the wave- 
length resolution of the stellar models. For every pixel in 
the output images, the attenuation, defined as the surface 
brightness in the image calculated including the effects 
of dust divided by the surface brightness in the image 
not including dust, is calculated for each of the radiative 
transfer wavelengths. This attenuation is then interpo- 
lated over all the wavelengths in the stellar-model SED, 
and this interpolated function is multiplied by the dust- 
free SED. This procedure yields an SED which has all 
the spectral features of the stellar model, with a large- 
scale behavior determined by the dust attenuation. The 
only exception to this is when the attenuation is ~S> 1, 
i.e. there is more light in the images with dust than in 
those without. This frequently happens when light is scat- 
tered by dust in regions that contain no emission, fn these 
cases, the attenuation is of questionable meaning, and the 
SED is interpolated directly from the radiative-transfer 
results. (Even in the cases where the attenuation is < 1, 
the interpolated SED is an approximation since it does 
not account for the small-wavelength spectral features of 
light scattered into the line of sight, only for those present 
in the direct light. When sunrise is updated to use the 
polychromatic algorithm for the galaxy simulations, this 
will no longer be an issue.) These data cubes are then 
integrated over a number of filter passbands to generate 
broadband images of the simulated object. 

Also during the postprocessing stage, the total lumi- 
nosity absorbed in every grid cell is calculated by interpo- 
lating the absorbed luminosity at the radiative-transfer 
wavelengths and integrating over all wavelengths. By en- 
ergy conservation this is equal to the total luminosity 
reradiated by the dust in the mid- and far-infrared, and 
images of the reradiated luminosity are then generated 
by assuming zero opacity and using Equation 1121 

While there is no self-consistent calculation of the 
dust temperature distribution and hence the SED of 
the radiation emit ted by the dust, as done by e.g. 
Missel t et alJ 1200 IT) , a rough idea of what the infrared 
SED should look like, in the context of galaxies , is pro- 
vided by the templates of iDevriendt et al.l 1119991) . These 
templates provide the dust emission SED of the galaxy 
as a function of its total infrared luminosity, which makes 
it possible to estimate the total luminosity in any given 
far-infrared passband. ft does not provide an estimate 
of the spatial variations of the dust-emission SED, so 
constructing images of the far-infrared emission can only 
be done assuming all cells have the same dust-emission 
SED. A self-consistent calculation of dust temperature is 
a planned upgrade to the code. 

Images are also generated for other quantities from 



5 TESTS 

With any scientific code, it is crucial to test it against 
independent results, either calculated analytically or us- 
ing another (hopefully correct) code. The full radiative- 
transfer problem, with nonisotropic scattering, absorp- 
tion and arbitrary geometry is far too complicated for 
analytic solutions, so the code has been tested against 
published results obtained with other Monte-Carlo codes. 
Sections 15.11 - 15.31 show comparisons with simpler sce- 
narios in which the geometry is well-defined and does 
not factor into the uncertainty. Section 15.41 show how 
the polychromatic algorithm performs. Finally, in Sec- 
tion E3| sunr ise is used to c a lculat e the 2D benchmark 
problem from IPascucci et al.l <I2004T) . which is also used 
to demonstrate the relative efficiencies of the mono- and 
polychromatic methods. 



5.1 Comparison to Witt (1977) 

In a series of papers, IWitd ll977allblcl hereafter collec- 
tively referred to as W77) used a Monte Carlo radiative 
transfer code to calculate the emerging surface-brightness 
distribution from a reflection nebula. This work is a suit- 
able test for radiative-transfer calculations, because the 
problem geometry is simple and the papers contain an ex- 
tensive study of how the surface-brightness distribution 
is affected by changing the input parameters. The nebula 
consists of a cylindrical, homogenous slab of dust. The 
star is located along the centerline of the cylinder, either 
in front of the dust, or embedded within it. The surface- 
brightness distribution was then described by the quan- 
tity S/F, where 5* is the surface brightness of the nebula 
and F is the flux from the star, observed at the distance 
of 126pc, as a function of the angular distance from the 
star. In sunrise, the adaptive grid was used to approx- 
imate the analytically defined cylinder used by W77. In 
practice, the results were insensitive to the degree of grid 
refinement used. 

A selection of the cases calculated by W77, varying 
the viewing angle, grain albedo, scattering asymmetry, 
and optical depth of the nebula, are presented in Fig- 
ure [5] In general, the results agree well, after a constant 
factor of 2 discrepancy in surface brightness has been ac- 
counted for. This discrepancy, in the sense that the W77 
results are too low, appears attributable to a normal- 
ization error in the W77 results (A. Witt 2000, private 
communication) , as it has been confirmed by other codes 
(A. Watson 2000, private communication). There is a sys- 
tematic trend where sunrise predicts a less steep rise in 
surface brightness close to the star, but this is consistent 
with the expected effect from the wide radial bins used 
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Figure 2. The radial brightness profiles of a nebula with an embedded star, in comparison with IWit J il977oh . The plots show 
the surface brightness of the nebula, normalized by the flux of the star, as a function of the angular distance away from the star. 
The SUNRISE results are shown as lines, while the W77c results (after multiplying them by 2, see text) are shown as points. The 
four panels explore the dependence of the results on four different parameters: In the upper left, the dependence on viewing angle 
(W77c, Table 7); in the upper right, dependence on scattering asymmetry (W77c, Table 10); in the lower left, dependence on 
albedo (W77c, Table 11); and in the lower right, dependence on optical depth (W77c, Table 4). In general, the agreement is good 
except that SUNRISE seems to systematically predict a less steep rise in brightness close to the star. As explained in the text, this 
is consistent with the wide bins used by W77 along with the steep rise in surface brightness for small radii. There also appears 
to be a trend, in the upper left plot, for the SUNRISE results to rise faster than the W77c data towards more oblique angles. The 
most seriously discrepant case is the a = 1 case in the lower left, where W77c predicts a much higher surface brightness than that 
predicted by SUNRISE. The origin of this discrepancy is unknown. 



by W77 in combination with the steep non-linear rise in 
surface brightness towards the center. 



The largest discrepancy is the a = 1.0 case in the 
lower left of Figure |5] where W77 predicts a surface 
brightness significantly higher than sunrise. The origin 
of this discrepancy is unknown, but unlike the W77 data, 
the sunrise results appear to be the extrapolation of the 
lower albedo data. 



5.2 Comparison to Watson & Henney (2001) 



IWatson fc Hennevl 11200 J) presented a problem geometry 
very similar to the W77 study, and argued that this was 
a suitable test for radiative-transfer codes in that it was 
simple but yet nontrivial. In this problem, an infinite 
plane-parallel slab of unit optical depth is illuminated 
either by a point source on the surface of the slab or 
by a collimated beam incident along the normal of the 
slab. (While the extent of the slab is described as infinite, 
this violates those authors' own definition of the set of 
problems they are treating, in which the opaque region is 
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Figure 3. Comparison with IWatson &: Hennevl feOOll) . showing the luminosity (or, more accurately, the radiative intensity) 
emerging in different directions from an infinite slab. On the left, the equivalent of their figure 1, where a point source is sitting 
on the surface of the slab, and on the right the equivalent of their figure 2, where a collimated beam is radiating into the surface. 
The WH01 results are shown as diamonds, the SUNRISE results as asterisks. Direct and singly scattered light is in blue (WH01 
separated this into direct and singly scattered light, but that is not easy to do by default in sunrise), double-scattered light in 
red and triply scattered and above in green. The results agree to a remarkable degree. 



only of finite extent. This restriction applies also to SUN- 
RISE, so the slab was made finite bu t much larger than 
its thi ckness.) In contrast with W77, IWatson fc Hennevl 
(2001) used the "radiative intensity" instead of the sur- 
face brightness. The radiative intensity is the total flux 
resulting from the source region, at large distance, mul- 
tiplied by the distance squared. This quantity, which is 
independent of distance, can be determined even in cases 
where the source is unresolved. Figure [3] shows the ra- 
diative intensity in different directions using the point 
source and collimated beam. The results have been split 
into three categories: direct plus single-scattered light, 
doubly scattered light, and more-than-doubly scattered 
light, (sunrise has no explicit facilities for selecting light 
depending on its history, but by comparing runs with 
different J m i n , this can be extracted.) The sunrise re- 
sults agree remarkably w ell with the ones presented by 
IWatson fc Hennevl fcOOllb 



5.3 Comparison to Witt & Gordon (1996) 

As a more comp licated test, the "clum py scattering 
environments" of IWitt fc Gordonl (1199a . from now on 
WG96), was replicated and the results compared. Briefly, 
WG96 investigated the escape of radiation from a clumpy 
environment consisting of a spherical region built up by 
smaller cubes of randomly assigned high- or low-density 
material illuminated by a central point source. The setup 
is described by 4 parameters: the filling factor, //, the 
fraction of cells which have high density; the density con- 
trast, fc, the ratio of opacities between the low-density 
and high-density cells; the grid resolution, N; and the 
optical depth of a homogeneous distribution with the 
same amount of dust, th- WG96 also defined the effec- 



tive optical depth t cS = — ln(I/ dir /L t ot), where L dir is 
the direct luminosity escaping through the distribution, 
without scattering, and L to t is the source luminosity. An 
image of a realization of these clumpy regions is shown 
in Figure 2] 

This test is more difficult to perform conclusively, 
because the problem geometry itself is random. To get 
accurate results, not only must the emerging luminos- 
ity be averaged over sufficiently many lines of sight, but 
enough realizations of the clumpy medium must also be 
used. WG96 did not described the number of lines of 
sight or the number of realizations of the medium used, 
but our results indicate that at least 150 lines of sight and 
60 realizations were necessary to obtain reasonably con- 
verged results. Even so, significant uncertainties in the 
SUNRlSEresults remain. The uncertainties in the results 
obtained by WG96 are unknown. 

Due to the nature of the grid used by sunrise, the 
problem cannot be compared exactly; in WG96, the ray 
tracing was terminated when the rays pass through a 
sphere inscribed within the grid, but this cannot be done 
in sunrise. Instead, the adaptive grid was used to ap- 
proximate a spherical surface, leading to a slightly "tiled" 
surface, as opposed to smooth, spherical one. Four sub- 
divisions were used, for which the results were largely 
converged. In any case, the uncertainty due to this effect 
is dwarfed by the uncertainty due to the averaging over 
lines of sight and realizations of the medium. 

Two tests are done. The first is the equivalent of fig- 
ures 8 and 9 in WG96, where the ratio of the scattered 
luminosity, Z/ S cat, the luminosity escaping the system af- 
ter scattering at least once, to the total luminosity is 
plotted against r e ft for varying values of the density con- 
trast fc. The results are shown in Figure |S] As the density 
contrast departs from unity, the effective optical depth 
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Figure 5. Comparison with figures 8-9 of lWitt &: Gordonl Jl996f) : scattered luminosity emerging from the clumpy system vs. 
effective optical depth, for a range of density ratios (1 > fea/fei > 0.001) in steps of factors of vlO, starting with k = 1 at r e g = 10. 
On the left, the equivalent of figure 8 in WG96, three different filling factors, //, are shown in different colors, while on the 
right, the equivalent of figure 9 in WG96, the colors indicate different grid grid resolutions, N. The WG96 points are connected 
by lines, while the corresponding SUNRISE points, for a number of different random realizations, are found immediately to the 
right and below. In general, as the density contrast departs from unity, the effective optical depth is lowered, and the scattered 
luminosity increases. The SUNRISE results agree with WG96 for the fc = 1 case, and for large N, but in other cases SUNRISE predicts 
systematically higher r e g and lower L sca t. 



is lowered and the scattered luminosity increases. This is 
a natural consequence of clumping, as low-density paths 
open up through the medium and allow radiation to es- 
cape. The sunrise results agree with those of WG96 very 
well for the homogenous case, but as the density contrast 
increases, the SUNRISE points show systematically higher 
T e flF than the WG96 ones. The discrepancy seems to de- 
crease as the number of grid cells increase; the TV = 40 
case agrees within the statistical uncertainty. 

The statistical spread in the sunrise points is 
strongly correlated. The main cause of this uncertainty 
is the random variation in the number of high- vs. low- 
density cells in the grid. Realizations in which the number 
of high-density cells is small will have smaller r e B and 
larger L sca t. Interestingly enough, this variation is also 
generally in the direction towards the WG96 results. Ex- 
periments showed that artificially lowering // by about 
0.015 from the // = 0.10 case, while keeping the opacities 
unchanged, will largely remove the discrepancy between 
the sunrise and WG96 points. However, this does not 
work very well for k < 0.01, and is probably not the cause 
of the discrepancy, the origin of which remains unknown. 

The second test is the equivalent of figure 15 in 
WG96. This compares the radial surface brightness dis- 
tribution of scattered light, averaged over many viewing 
directions and realizations of the clumpy medium, for the 
case tb = 10, // = 0.10 and N = 10. The results are 
shown in Figure HJ As the density contrast departs from 
unity, the overall surface brightness of the scattered light 
increases from the homogenous case, and the distribu- 
tion becomes more peaked towards the center. For very 
small contrast ratios, when the low-density cells are al- 



ready optically thin, the surface brightness again begins 
to decrease. 

Because the WG96 figure is plotted with "arbitrary 
units" on the x-axis, it has been assumed that the inter- 
val plotted is the full extent of the object. The units of 
their surface brightness are also not specified, so the SUN- 
RISE results have been adjusted so that the k = 1 case 
agrees with WG96. The sunrise lines agree well with 
the WG96 results over a large range of radii, but fall off 
more quickly towards the edge. The sunrise results are 
also more sharply peaked compared to the innermost bin 
of WG96, but this is because the direct light from the 
point source in the center has not been removed. There 
is also agreement in some of the finer structure, for ex- 
ample a small bump in the surface brightness at r ~ 450 
for the k — 0.001 case, corresponding to the radius of the 
innermost grid cell boundary. 

As seen, there are some systematic differences be- 
tween the WG96 results and those determined using 
SUNRISE. However, the discrepancies are not large, and 
both results share the same qualitative behavior. While 
it would be desirable to pin down the cause of the dis- 
crepancies, we do not believe that they are severe enough 
to cause alarm. 



5.4 Tests of the Polychromatic algorithm 

The "polychromatic" algorithm has been tested with a 
preliminary implementation, showing that it does indeed 
reproduce the results obtained with the single- wavelength 
calculation, but it is not yet included in the production 
version of sunrise. 

As a test, the W77 comparisons presented in Sec- 
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Figure 7. The polychromatic algorithm tested on the W77 problems. The three plots show the deviation in surface brightness 
between the polychromatic results, where all 13 test cases were calculated simultaneously, and the monochromatic results presented 
in Figure [5] The different cases have been displaced from zero for clarity. For each test case, results using three different reference 
opacities and scattering asymmetries are shown by different colored lines. In general, the agreement is good. There is significantly 
more noise in those cases where the reference values differ greatly from the case being tested, most notably the bottom blue line in 
the left plot where very forward-scattering dust (g = 0.9) is being traced with practically isotropically scattering dust (g rc f = 0.1). 
Increased noise is also evident in the top left green line (g = 0.1, g IC { = 0.9) and the bottom right blue line (k = 5, K rc f = 1). 
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Figure 8. Spectra from the galaxy merger simulations, exemplifying the effects of the new polychromatic algorithm. In both 
images, the red line shows the intrinsic emission in the pixel, neglecting radiative-transfer effects. The green line shows the 
old algorithm, where the spectrum is interpolated between 20 discrete wavelengths, and the blue line is the result from the 
polychromatic radiative transfer. The top plot shows a pixel containing an obscured H II region, the bottom plot a pixel near 
an Hll region where the UV flux is dominated by scattered light. While the results agree well on the overall spectral shape, the 
new method gives significantly more realistic results especially for the small-scale spectral features, at a fraction of the runtime. 
Note in particular how the polychromatic algorithm predicts the appearance of the Ly a absorption line, the disappearance of the 
nebular Balmer continuum edge, and the increased attenuation of the Paschen/3 line at 1.3 /im in the spectrum of the Hll region. 
The polychromatic calculation, including 500 wavelengths, used about 8 times the CPU time required for one monochromatic 
calculation. With monochromatic ray tracing, 500 separate wavelengths would have to be used to predict the same amount of 
detail, which would require a factor of 50 more CPU time. 
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Figure 4. Image ol the radiation emerging from a realization 
ol the Witt & Gordon (1996) random clumpy medium, with 
parameters t h = 10, // = 0.10, k = 0.01, and N = 20. 
Most radiation is emerging on paths through the low-density 
medium. The inner surfaces of the high-density clumps can be 
seen as bright edges due to the efficient scattering of light off 
of them, while dark shadows are cast towards the outside. 



tion l5.1l were redone, but now all different optical depths, 
albedos and scattering asymmetries were calculated si- 
multaneously. The results are shown in Figure On 
the whole, the agreement is excellent when the reference 
opacity and scattering asymmetry is close to the case be- 
ing calculated. When the reference parameters are very 
different from the case being calculated, the results can 
indeed become very noisy. This is most evident in the 
cases where very forward-scattering dust (g = 0.9) is be- 
ing calculated using an isotropically-scattering reference 
case (g rc f — 0.1, or vice versa. In these cases, the weight- 
ing factor can reach values of up to 200. 

Very large weighting factors can also be obtained in 
very optically thick cases where t /r rc f < 1. In fact, it is 
clear from eauation l^ll that in those cases w\ will increase 
without bound while the total contribution w\exp{— r rB f) 
from those cases remains finite. This suggests that con- 
vergence will be poor, and emphasizes that a proper 
choice of reference parameters is crucial for the accuracy 
of the method; the reference wavelength should be chosen 
such that the range of weighting factors encountered in 
the problem is minimized. 

The polychromatic algorithm has also been tested on 
our galaxy merger simulations. This is a favorable situ- 
ation, as direct light often dominates. Examples of out- 
puts are shown in Figure [S] which shows the spectra of 
two different pixels in a galaxy merger simulations. 

The dust opacity in a standard Milky- Way dust 
model varies by close to three orders of magnitude from 
the far-UV to the near-IR. The disadvantages of the 
large weighting factors resulting from this large range of 
opacity can be alleviated by stratifying the calculation 
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Figure 6. Radial surface brightness distributions of light 
for a range of density ratios (1 > &2/fci > 0.001) . This is 
the equivalent of figure 15 in W itt fc Gordon! (1996J), but for 
clarity the different cases have been separated into two plots. 
The SUNRISE results are shown as lines, and the WG96 results 
as points. The SUNRISE results fall off more sharply towards the 
edge, but over a large range of radii the agreement is good. 
The sharp rise in the SUNRISE lines compared to the innermost 
WG96 point is because the direct light from the central point 
source is included in the SUNRISE results, but has been removed 
in the WG96 data. 



into several ranges of wavelengths, which will restrict the 
range of opacities treated in the same random walk. This 
is used in the next section. 



5.5 The Pascucci et al. (2004) 2D benchmark 

IPascucci et alJ B2004L from now on P04) aimed at estab- 
lishing a standard benchmark radiative-transfer problem, 
in the spirit of similar benchmark hydrodynamic tests, 
where a number of codes are compared against each other 
in a more complicated problem that lacks an analyti- 
cal solution. They designed a two-dimensional, axisym- 
metric problem modeling a circumstellar disk and cal- 



© 2006 RAS, MNRAS OOO.IT1I231 



18 Patrik Jonsson 



,=0.1 



s? 



0.02 
0.01 
0.00 



01 
02 

02 
01 
00 
01 
02 

02 
01 
00 
01 
02 

10- 



Polychromatic 

Monochromatic 


1=12.5°- 


- 


- 


■ 


i=42.5°- 


- 


- 


: — - A ^^______ 


i = 77.5°- 


- 


- 



0.2 

0.0 
-0.2 

0.2 

0.0 
-0.2 

0.2 

0.0 

-0.2 



io- 



10- 
A/m 

T„= 1 



Polychromatic 
Monochromatic 



i = 1 2.5° 



i = 42. 




10- 



fr? 



11 







T v = 






0.1 
0.0 
0.1 




Polychromatic 

Monochromatic 


i = 12.5 


-i 








1 


0.1 
0.0 

0.1 






i = 42.5 


i 








i 


0.1 
0.0 
0.1 






i = 77.5 


-I 








-I 



10- 



10- 
A/m 

T v =100 



Polychromatic 
Monochromatic 



^12.5° 



1=42. 



Stratified 



i = 77.5° 




A/m 



10- 



A/m 



Figur e 9. Difference in SED between the results obtained with SUNRISE and the reference code RADICAL for the lPascucci et al.l 
C2004) benchmark problem. The monochromatic algorithm is shown in red, the polychromatic algorithm in blue. For the r v = 100 
edge-on case, the "stratified" polychromatic calculation (see text) is shown in green. The error bars (frequently too small to be 
visible) indicate the lcr Monte-Carlo sampling uncertainty in runs with 10 6 rays (per wavelength for the monochromatic algorithm). 
The discrepancy reaches 40% for the edge-on r = 100 case. The results of the polychromatic algorithm are indistinguishable from 
those obtained with the monochromatic one except for the most optically thick case, where the results diverge around 1 /im. 
The stratified polychroma tic calcu lation avoids this problem and agrees well with the monochromatic run. This figure is directly 
comparable to Figure 8 in lPascucci et alJ J2004). 



culated the dust temperature distribution and emerging 
SED from a number of inclination angles. The optical 
data for the dust grains and the outputs from the 5 codes 
they used are available for download, which makes a de- 
tailed comparison possible. Since sunrise currently does 
not have the capability to calculate temperature distribu- 
tions, only the SEDs at wavelengths shorter than 5/xm, 
where stellar light dominates, were compared. Compar- 
isons were done both using the monochromatic and poly- 
chromatic algorithms, which makes it possible to evalu- 
ate the efficiency of the polychromatic method in a fairly 
complicated problem. 

The codes tested by P04 all used two-dimensional, 
axisymmetric grids with variable grid spacings. The 



orthogonal, explicitly three-dimensional adaptive-mesh 
grid used by sunrise is not ideally suited for such a 
problem, as the number of grid cells necessary to re- 
solve the problem geometry well is much greater com- 
pared to a 2D grid. T he P04 benchmark was also used by 
lErcolano et al. (2005) to test the MOCCASIN combined 
photoionization/dust radiative transfer code, which uses 
a (uniform) Cartesian grid. 

The sunrise results are presented in Figure^] which 
is analogous to Figure 8 in P04. The differences be- 
tween the results using su nrise and the RADICAL 
iJDullemond fc Turoll a 2000) code are plotted as a func- 
tion of wavelength for the four different optical depths 
and three different inclination angles used by P04. For 
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Figure 10. The efficiencies of the polychroma tic and 
monochromatic methods in the IPascucci et alJ fe004l) bench- 
mark problem. The polychromatic algorithm is shown as solid 
lines, the monochromatic as dashed lines, while the color in- 
dicates the optical depth of the problem. For the r v = 100 
case, the efficiency of the "stratified" polychromatic calcula- 
tion is shown as a dot-dashed line. For lower optical depths, 
the efficiency of the polychromatic algorithm exceeds that of 
the monochromatic one for all wavelengths. For t v = 10, the 
minimum efficiency, indicating the error of the most poorly 
constrained wavelength, is still higher for the polychromatic 
algorithm. For the highest optical depth, the polychromatic 
calculation suffers from the large range of optical depths for 
different wavelengths, and shows very low efficiencies around 
0.8 /im. The stratified calculation, where wavelengths longer 
than 0.5 (im are traced separately, avoids this problem and 
keeps the minimum efficiency significantly higher than the 
monochromatic algorithm. 



small optical depths the results agree very well, which 
is not surprising since the SED is close to the intrinsic 
blackbody SED of the central source. For larger optical 
depths, and especially for the edge-on configurations, the 
discrepancies reach ±40%. This is larger than the internal 
differences between the codes used by P04, and is prob- 
ably due to the use of an orthogonal three-dimensional 
grid. This hypothesis is supported by the fact that the 
results were converging (slowly) towards the P04 results 
as the grid resolution was increased. The results plot- 
ted used a grid with 1.7 x 10 6 cells and a minimum 
cell size of 0.02 x 0.02 x 0.002 AU . Grids with up to 
4.3 x 10 6 cells have been tried, and improved the agree- 
ment of around 5% in the most optically thick case. The 
agreement between the P04 benchmark and the sunrise 
res ults is still sig nificantly better than those presented 
bv lErcolano et all (2005ft for the MOCCASIN code, us- 
ing a uniform Cartesian grid, which disagree with the 
P04 results by up to a fac tor of 20. It is unclea r what 
grid resolution was used bv lErcolano et al. (2005). 

The benchmark was calculated using both the 
monochromatic and polychromatic methods. For the 
monochromatic calculations, 10 6 rays were used for each 
of the 28 wavelengths used by P04 blueward of 5 pm, 



while the polychromatic calculations used 10 6 polychro- 
matic rays and a reference wavelength of 0.4 pm. Except 
for the r v = 100 case, the results are indistinguishable 
in Figure As the optical depth increases, the poly- 
chromatic calculation becomes increasingly affected by 
the large range of optical depths encountered at different 
wavelengths. The reference wavelength used was 0.4 pm, 
where both opacity and albedo are high, and as a con- 
sequence accuracy suffers at wavelengths between 0.7 pm 
and 1.5 /im, where the opacity drops but the albedo is 
still fairly high. This problem was solved by "stratifying" 
the calculation into two wavelength ranges. The reference 
wavelength was kept at 0.4 pm for wavelengths shorter 
than 0.55 pm, while longer wavelengths were traced sep- 
arately, using a reference wavelength of 0.7 pm. This two- 
step calculation agrees very well with the monochromatic 
results. 

To quantify the relative efficiencies of the two meth- 
ods, the "efficiency", e, defined as 



(39) 



K 



where T is the CPU time required to complete the calcu- 
lation and <jf x is the Monte Carlo sampling uncertainty 
in the SED, is calculated. The efficiency defined in this 
way is insensitive to the number of rays traced and quan- 
tifies the (inverse of the) CPU time necessary to produce 
results of unit relative accuracy. Figure UTil plots the effi- 
ciencies of the monochromatic and polychromatic meth- 
ods as a function of wavelength for the edge-on cases 
in Figure |U] No effort was made to optimize the relative 
number of rays at different wavelengths in the monochro- 
matic calculation. 

For t v < 10, the efficiency of the polychromatic 
method exceeds that of the monochromatic one for all 
wavelengths. As mentioned above, the polychromatic al- 
gorithm begins to suffer at higher optical depths. For 
t v — 10, the efficiency of the monochromatic algorithm 
overtakes the polychromatic one at wavelengths longer 
than 0.5 pm. However, the minimum efficiency, which 
corresponds to the wavelength with the largest error in 
the SED, is still higher for the polychromatic algorithm. 

For the very optically thick r v = 100 case, the ef- 
ficiency of the polychromatic algorithm plummets for 
the same wavelength range where the errors were large 
in Figure [5] reaching a minimum value much lower 
than the monochromatic calculation. The stratified poly- 
chromatic calculation, however, maintains an efficiency 
greater than the monochromatic calculation out to 1 pm 
and has a minimum efficiency significantly larger than 
the monochromatic results. 

It is important to remember that this example used 
a comparably low number of 28 wavelengths. While the 
run time for the monochromatic method will increase lin- 
early as the wavelength resolution is improved, the poly- 
chromatic method will scale much better (at least until 
the run time becomes dominated by the vector calcula- 
tions as opposed to the ray tracing). It is clear that the 
particular strength of the polychromatic algorithm is in 
processing densely sampled wavelengths, not in covering 
a large range in wavelength. 
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Figure 11. The memory access pattern in a test run of SUN- 
RISE. The graph shows the cumulative distribution of strides 
(differences in memory address) when accessing consecutive 
grid cells during a run consisting of about 10 5 grid-cell-to-grid- 
cell steps. The minimum stride is 48 bytes, the size of the grid- 
ccll data structure, and this is the source of the discretization 
of the small strides. The vertical lines indicate 4096 and 16384 
bytes, the size of the TLB entries on the POWER3/Opteron 
and Altix 300 processors, respectively. The Hilbert ordering re- 
sults in noticeably smaller strides than the default depth-first 
ordering, but not enough to outweigh its increased complexity. 



6 CODE IMPLEMENTATION DETAILS 

SUNRISE is written in C++, with a highly modular code 
structure. It is easy to add new types of emission, sources 
of scattering, etc., as these are isolated from the core 
radiative-transfer engine. Efficient vector calculatio ns are 
provided by the Blitz++ library JVeldhuizenlll99ct) . 

SUNRISE uses the implementation of the 
Mersen ne Twister MT19937 rando m- number algo- 
rithm JMatsum oto fc Nishimura|l l99cfl provided by the 
Blitz++ library JVeldhuizenlll99ck ~ 



6.1 Grid Traversal 

The most time-consuming part of the calculation is the 
grid traversal, walking the ray from grid cell to grid cell. 
Because the rays can traverse the grid in any direction, 
the memory access pattern of the code is hard to predict, 
and often not optimal. To minimize CPU cache and TLB 
(Translation Lookaside Buffer) misses, it is important to 
access memory in as localized a fashion as possible. This 
is easy in applications like matrix multiplications, where 
it is up to the programmer to decide in which order to 
access memory, but hard in the application considered 
here. The solution is to order the grid cells in memory 
in such a way that locality is preserved; grid cells that 
are close in space should also be stored close in memory. 
Since this is a mapping from 3D to ID, it cannot be ac- 
complished perfectly, and various approaches exist. The 
strategy used by sunrise is to simply store the subgrids 
consecutively in memory in the order they are encoun- 



tered when traversing the tree of refined cells in a depth- 
first fashion. This actually results in a fairly compact 
memory access pattern. To explore the potential for im- 
provements, storing the grid cells according to a Hilbert 
curve was also tried, using t he algorithm described in 
iBartholdi fc Goldsmanl (1200 ll) . The Hilbert curve is an 
example of what is known as "space- filling curves" , a one- 
dimensional curve which fills three-dimensional space, 
and has b een argued as being nearly op timal for these 
problems fNied ermeier fc Sand ers 1996). The result of 
this comparison is shown in Figure fTTl While the Hilbert 
ordering did have a more compact access pattern, the dif- 
ferences were comparatively small and did not outweigh 
the accompanied increased complexity. Thus, the default 
memory ordering was retained. 

In addition, finding neighboring cells in the grid is 
quite expensive, as the tree has to be traversed up the hi- 
erarchy and then down again, accessing many higher-level 
grid cells in the process. Because of this, sunrise utilizes 
a caching scheme whereby each cell saves pointers to its 
neighboring cells once they have been determined. This 
speeds up subsequent neighbor searches considerably. 



6.2 Runtime Requirements 

The most obvious scaling is that the runtime of the 
ray tracing (excluding initialization, I/O, etc.) obviously 
scales linearly with the number of rays traced. The time 
required per ray is affected by many factors, including the 
problem geometry, the importance of scattering, and the 
terminating intensity I m i n . The time is dominated by the 
tracing of the ray from one grid cell to the other, which 
takes approximately constant time. Emission and inter- 
action events are much more rare and do not constitute 
the dominant computational load. From this it would be 
expected that the time required to trace one ray would 
scale with the number of grid cells, which is the case. 

1/3 

The scaling is roughly i r ~ N c , which would also be 
the naive expectation since this is how the number of 
grid cells traversed by a ray should scale with the total 
number of cells. 

Furthermore, increasing the number of cameras in- 
creases the computational load, since part of the ray 
propagation consists of calculating r obs , which means 
traversing the cells from the point of emission or interac- 
tion to each of the observers. The time taken to complete 
one Monte-Carlo ray can thus be split into two compo- 
nents such that t r — to + N c t c , where to is the time taken 
to complete the ray tracing in the absence of observers, 
and t c the time taken to traced the optical depth to an 
observer. This relation is indeed obeyed by the code, and 
the two components to and t c are such that, in typical 
cases of simulated galaxies, the time proportional to N c 
becomes dominant for N c > 3. 



6.3 Parallelization 

Any high-performance numerical code relies on parallel 
execution to reach high performance, so this capability 
was a basic requirement for the development of SUN- 
RISE. The Monte-Carlo method is easy to parallelize; as 
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long as the entire grid fits in memory, each ray is inde- 
pendent of others, sunrise uses multithreaded (shared- 
memory) parallelism during the grid creation and ray- 
tracing stages. As every ray is independent, communi- 
cation requirements are minimal resulting in essentially 
perfect scaling as long as the memory bus is not satu- 
rated. 

Completing the full calculation for one simulation 
snapshot takes 6-12 CPU hours on a contemporary dual- 
CPU 2.2GHz Opteron system, for typical conditions of 
11 camera positions, 21 wavelengths and 10 6 rays per 
wavelength. The total amount of CPU time consumed 
by the radiative-transfer calculations is roughly equal to 
the time used by the hydrodynamic galaxy merger simu- 
lations. 

On the NASA Altix 3000 system Columbia, a cluster 
of machines each consisting of 256 two-processor nodes, 
where the full, system-wide memory is visible as a sin- 
gle address space using high-speed interconnects, sunrise 
can be used on up to 16 processors with a penalty of only 
23 percent compared to running on one node. 

Adding the ability to run using distributed memory 
would vastly increase the size of the problems that could 
be treated. However, it would also make the paralleliza- 
tion much more complicated, and as the current ability is 
sufficient for our computational needs, such an upgrade 
is not currently planned. 



6.4 Distribution 

As a service to the community, sunrise is available 
under the terms of the GNU General Public License 
llFree Software Foundation! I1991I) 1 . Other users are en- 
couraged to use SUNRISE for their radiative-transfer ap- 
plications, and to add enhancements that increase its ca- 
pabilities. 



7 FUTURE IMPROVEMENTS 

There are numerous improvements that can be made to 
the radiative-transfer code. The polychromatic algorithm 
is currently being implemented in the production version 
of sunrise, and schemes for minimizing the impact of 
large weighting factors will be explored. One possibility 
is to split rays wit h large weights and sca tter them in dif- 
ferent directions IIDupree fc Fr alcv 2002J), which will in- 
crease the sampling in the heavily weighted part of phase 
space. Such a scheme could be combined with "Russian 
roulette" into a scheme where the code attempts to keep 
the weights of all rays within some specified range. 

The biasing of path lengths done in the polychro- 
matic algorithm opens up other possibilities; it is well 
known that the Monte Carlo method has problems treat- 
ing very large optical depths, for example in the case of 
a clou d heated by external radiation studied by I.Tuvelal 
(2005), because few rays make it into the opaque regions. 



1 The SUNRISE source code, documentation, and exam- 
ple outputs are available on the SUNRISE web site at 
http : //sunrise . f amilj enj onsson . org 



ljuvelal J2005TI explored the effects of biasing the distribu- 
tions of the external photons and the scattering direction. 
Another possibility would be to bias the path lengths 
between scatterings to larger values, which would allow 
better sampling of the rare photons which penetrates un- 
usually deep into the cloud. 

These possibilities suggest that it would be advanta- 
geous to make it possible to not only use arbitrary scat- 
tering phase functions or emission distributions in SUN- 
RISE, but also to provide an infrastructure for customized 
"biasing plug-ins" which can be selected by the user de- 
pending on the problem at hand. 

Another obvious improvement is to include a self- 
consistent calculation of the dust emission, so predictions 
of infrared observations of galaxies with e.g. the Spitzer 
Space Telescope can be made. These calculations can be 
done to varying degrees of sophistication, from calculat- 
ing equilibrium temperatures of single grain species to full 
stochastic temperature distributions of very small grains 
and PAH emission, but it becomes necessary to include 
the heating of dust grains due to the emission from other 
grains. This introduces a coupling between the local emis- 
sivity and radiation intensity, which makes it necessary 
to integrate until the dust temperature distribution con- 
verges. 

An interesting alternative, applicable in situations 
where the opacity is not a function o f temperature, is the 
method of lBiorkman fc W ood (2001). In this method, the 
dust temperatures are updated whenever an absorption 
event happens and the energy is reemitted with an SED 
equal to the difference between the current emission SED 
and the one before the absorption event happened. This 
ensures that the net radiation emitted by the dust is ap- 
propriate for the temperature of the dust grains, and as 
more rays are traced the dust temperature distribution 
will relax toward t he equilibrium value without iteration. 
IfBaes et alJ IJ2005I) studied this method and concluded 
that while it produces the correct frequency distribution 
in the case of grains in thermal equilibrium, the method 
will not work for stochastically heated grains because the 
probability distribution required for the reemitted radi- 
ation will become negative. As negative probabilities are 
unphysical, the method will fail. The polychromatic algo- 
rithm suggests a way around this obstacle: There is no ob- 
vious reason why the weights of the rays in the polychro- 
matic algorithm can't be negative. The problem is not 
fundamentally that negative probabilities are necessary, 
but that the grain has emitted too much energy in a cer- 
tain part of the spectrum and this must be corrected for. 
This can be accomplished with negative weights for cer- 
tain wavelengths. Since it will only serve to remove energy 
which was previously emitted, the radiation field should 
still converge towards the true value, only it will not con- 
verge from below as in the original iBiorkman fc Wood! 
(2001) formulation. This possibility will be explored in 
SUNRISE in the future. 

Another important improvement for correctly pre- 
dicting the SED of galaxies is a more detailed modeling of 
star-forming regions, which are not resolved in our galaxy 
simulations. Studies have shown that extra attenuation 
of young stellar populati ons are necessary for fitting dust 
attenuations in galaxies JSilva et alll998ilCharlot fc FalJ 
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l200d: iTuffs et al]l2004h . While the adaptive grid could 
be used to resolve these regions, the computational cost 
would be prohibitive. A better way would be to use sub- 
resolution m£d^ls_of_tteemission from star-forming re- 
gions fe.g. lOroves et al.ll2004l and feed this emission into 
the sunrise grid. 



8 CONCLUSION 

This paper has described sunrise, a new Monte-Carlo 
code for calculating the radiative transfer of light through 
a scattering and absorbing medium, su nrise builds on 
advanced Monte-Carlo codes ( Kurosawa ct al. 



previous 



20041: [Baes et afllgOO.I: lOordon et alJl200lHBianchi et al 
200d:lLucvlll999l:IWolf et al.lll999l) . and adds a polychro- 
matic algorithm, where all wavelengths are traced simul- 
taneously, and efficient parallel computation in a flexi- 
ble, modular package, making calculations with a spatial 
dynamical range of more than 10 4 feasible. Images at 
any wavelength from far-UV to near-IR from an arbi- 
trary number of directions are generated, as well as the 
radiation intensity and dust luminosity as a function of 
position in the object. In addition, sunrise includes a 
framework for calculating the effects of dust in hydrody- 
namic simulations, and the code is freely available to the 
community. 

Accurate radiative-transfer calculations with realis- 
tic geometries are crucial in tying theory to observations 
in any situation where dust significantly affects the ra- 
diation emerging from an object. This applies to such 
diverse a family of objects as galaxies, star-forming re- 
gions, AGN, and protostellar disks. The outputs from the 
radiative-transfer calculations generate "simulated obser- 
vations" of the object, so comparisons with observations 
can be done using observational instead of theoretical 
quantities. This approach requires specifying many free 
parameters, but in reality only serves to make this depen- 
dence explicit, since if the comparisons are done with the- 
oretical quantities, the same parameters normally need 
to be specified to convert observed quantities to intrinsic 
ones. In many cases, such as the complicated radiative- 
transfer situations that can be solved by sunrise, it is 
not even possible to invert the observations, and making 
observational predictions from a theoretical model is then 
the only avenue possible. Another advantage of compar- 
ing observational quantities is that it is usually easier to 
mimic selection biases and instrumental effects in simula- 
tions than to infer their effects on theoretical quantities. 

sunrise is currently being used to investigate dust 
attenuation in simulations of isolated spiral galaxies 
(Rocha Gaso et al., in preparation), to generate a library 
of images of simulated galaxy mergers which can be di- 
rectly compared to Hubble Space Telescope observations, 
and to quantify the timescales over which mergers result 
in disturbed morphologies. 
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